function [P Pi c] = AIS_fixedpoint_hat(Y,Xi,that,t)

% Ouput: Price indices P and Pi
% Pi is Pi_hat^-theta, P is P_hat^-theta

N = length(Y);
YW = sum(Y);
Pi = ones(N,1);
P = ones(N,1);
lim = 1e-8;
c = 0; diff=1;
while diff>lim

  Pinew = Y./(Y+sum(Xi,2)).*Pi.^(1/(t+1))./P + sum(Xi./repmat(Y+sum(Xi,2),1,N).*that.^(-t).*repmat(Pi',N,1),2);
  
  % Rescale
  Yhatnew = Pinew.^(1/(t+1));
  Pinew = Pinew/sum((Y/YW).*Yhatnew).^(t+1);       %     
  
  Pnew = Y./(Y+sum(Xi,1)').*Pinew.^(1/(t+1))./Pinew + sum(Xi./repmat(Y'+sum(Xi,1),N,1).*that.^(-t).*repmat(P,1,N),1)';    

  x = Pnew-P;  
  diff = x'*x;
  Pi = Pinew; 
  P = Pnew;
  c = c+1;
  if c>10000 break; end;
end
disp('Converged. Counter:'); disp(c); disp('Latest diff'); disp(diff);
